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1 Introduction 


Predictive capabilities of numerical tools for progressive damage and failure analy- 
sis (PDFA) of composite structures are predicated upon the robustness, accuracy, 
and objectivity of the tools. As sophisticated numerical tools become more feasible 
for widespread use in industry, significant weight and cost saving in the design of 
lightweight composite structures will be apparent. Virtual testing of materials, incor- 
porating PDFA, can be used to evaluate the viability of materials or configurations 
prior to further scrutiny via physical testing. This grants analysts more flexibility dur- 
ing preliminary structural design stages and will ultimately manifest as more efficient 
and cost effective designs. 

Continuum damage mechanics (CDM) has emerged as a viable option for predict- 
ing the non-linear behavior of composite structures. The first CDM theory was de- 
veloped by Kachanov (1958, 1986). Subsequently, many publications on this subject 
were produced, including numerous books [ Talreja (1985a); Lemaitre and Chaboche 
(1994); Lemaitre (1996); Krajcinovic (1996); Voyiadjis and Kattan (2005)]. Typi- 
cally, a set of scalar damage variables, or internal state variables (ISVs), introduce 
anisotropic damage into the composite constituent behavior by penalizing the com- 
ponents of the material stiffness tensor, and non-linear functions are used to control 
the damage evolution. Various authors have used crack density, geometry, strain 
energy release rate, and other crack features to characterize the damage evolution 
[. Dvorak et al. (1985); Talreja (1985b); Laws and Dvorak (1988); Lee et al. (1989); 
Nairn (1989); Tan and Nuismer (1989); Gudmundson and Ostlund (1992); McCart- 
ney (1992a, 1998)]. Others postulate damage evolution laws and characterize those 
laws using experiments [Allen et al. (1987a, b); Talreja (1994); Paas et al. (1992); 
Matzenmiller et al. (1995); Bednarcyk et al. (2010)]. CDM models must also employ 
failure criteria to indicate damage initiation. More recently, increasingly sophisti- 
cated failure criteria have been developed to better represent the phenomenological 
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behavior of a damaging composite lamina [ Puck and Schilrmann (1998, 2002); Pinho 
et al. (2005)] and used in conjunction with CDM. 

CDM techniques offer computationally efficient and readily implementable means 
to capturing the effects of damage and failure in composite materials. Unfortunately, 
the majority of the criteria and evolution laws are formulated upon phenomenological 
observations of the directional dependence of damage evolution, rather than modeling 
the physics of the actual damage mechanisms. Separate damage variables are used to 
degrade different components of the stiffness tensor (directions) depending on whether 
the damage is said to accrue in the matrix or fiber constituents of the composite, but 
the variables do not explicitly distinguish between the separate damage mechanisms. 
Furthermore, many theories involve a multitude of parameters that are difficult to 
measure and must be calibrated to correlate with experimental data. 

When implemented within the finite element method (FEM), many PDFA method- 
ologies that utilize CDM breakdown when the material enters the post-peak strain 
softening regime locally within an element. Loss of positive definiteness of the tangent 
stiffness tensor leads to pathological mesh dependence [ Bazant and Cedolin (1979); 
Pietruszczak and Mroz (1981)]. To overcome this deficiency, Bazant (1982) devel- 
oped the smeared crack, or crack band, model that introduces a characteristic ele- 
ment length into the formulation of the damage evolution. The original formulation 
assumed that the mode I crack band always aligns with the principle axes [ Bazant 
and Oh (1983)]. de Borst and Nauta (1985) altered the formulation to accommo- 
date a fixed crack band under mixed mode conditions. An encompassing overview of 
smeared crack band models is provided by Spencer (2002). 

A thermodynamically-based, work potential theory, known as Schapery theory 
(ST), was developed for modeling matrix microdamage in fiber-reinforced laminates 
(FRLs) [ Lamborn and Schapery (1988); Schapery (1989, 1990); Lamborn and Schapery 
(1993)]. Sicking (1992) and Schapery and Sicking (1995) extended the formulation 
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to include the effects of transverse cracking by adding an additional internal state 
variable (ISV) and predicted the evolution of microdamage and transverse cracking 
in coupon laminates analytically. Pineda et al. (2010) implemented this extended 
formulation in a numerical setting to simulate the failure of a buffer strip-reinforced, 
center-notched panel (CNP). However, due to the cumbersome nature of the evolu- 
tion equations, the microdamage and transverse cracking evolution equations were 
decoupled to arrive at a more efficient implementation. Since no characteristic length 
is introduced into the formulation, the theory produces mesh-dependent results in a 
computational setting. 

The ST formulation is modified here, resulting in the enhanced Schapery theory 
(EST), to include the effects of macroscopic transverse and shear matrix cracking, as 
well as fiber breakage, using an approach that differs from Sicking (1992); Schapery 
and Sicking (1995); Pineda et al. (2010). A deliberate distinction between damage and 
failure is made. Damage is defined as the effects of any structural changes resulting in 
a non-linear response that preserves the positive definiteness of the tangent stiffness 
tensor of the material. Conversely, failure is considered to be the consequence of 
structural changes that cause post-peak strain softening in the stress versus strain 
response of the material. Here, matrix microdamage is categorized as a damage 
mechanism, but macroscopic matrix cracking and fiber breakage are hypothesized to 
be failure mechanisms resulting from damage localization. The traditional ISV used in 
ST is maintained to model microdamage. Upon failure initiation, the element domain 
is no longer considered a continuum, and a smeared crack approach is used to model 
the embedded discontinuities. Three new IS Vs, which incorporate the characteristic 
length of the finite element, dictate the evolution of the failure mechanisms. The EST 
formulation presented in Section 2 offers non-linear progressive damage coupled with 
mesh objective, post-peak strain softening. 

Mesh objectivity is demonstrated in Section 3. In Section 4, EST is verified against 
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experimental results for two center- notched panels (CNPs). Global load versus de- 
flection data, local strain gage data, as well as observed failure mechanisms obtained 
from experiments performed at the NASA Langley Research Center (LaRC) and ex- 
hibited in Bogert et al. (2006); Satyanarayana et al. (2007) are compared to numerical 
results. 


2 Enhanced Schapery Theory 

The previously developed ST ([ Schapery (1990, 1995); Schapery and Sicking (1995); 
Basu et al. (2006); Pineda et al. (2009, 2010)]) is extended to accommodate mesh 
objective, post peak strain softening. Separate IS Vs are used to govern the evolution 
of matrix microdamage, transverse (mode I) matrix failure, shear (mode II) matrix 
failure, and fiber breakage (mode I). The first and second laws of thermodynamics 
are enforced, establishing thermodynamically consistent evolution laws for progressive 
matrix microdamage, as well as post-peak failure. The following sections detail the 
formulation of this work potential theory. 

2.1 Thermodynamically-based Work Potential Framework 

As a material is loaded, a measure of the work potential facilitates modeling structural 
changes in the material, such as microcracking, which affect the elastic properties of 
the material. Energy that is not dissipated is recovered when the structure is un- 
loaded, and the magnitude of energy recovered is contingent upon the degraded, 
elastic properties at the previously attained maximum strain state. It is assumed, 
upon subsequent reloading, that the material behaves linearly, exhibiting the elas- 
tic properties observed during unloading, until the material reaches the preceding 
maximum strain state. After this state is achieved, structural changes resume, af- 
fecting degradation of the instantaneous elastic moduli of the material. This process 
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is shown in the uniaxial stress-strain curve displayed in Figure 1. The shaded area 
above the unloading line represents total dissipated potential Ws, and the triangular 
area underneath is the total elastic strain energy density We- It is assumed that the 
material behaves as a secant material and there is no permanent deformation upon 
unloading. This a reasonable assumption for FRLs Sicking (1992); however, plastic 
deformation can also be incorporated, if necessary [ Schapery (1990)]. Extension to 
treat viscoelastic and viscoplastic response is outlined in Hinterhoelzl and Schapery 


(2004). 

Both We and Ws are functions of a set of ISVs, S mi (m — 1,2 , M). These ISVs 
account for any inelastic structural changes in the material. Differentiating Ws with 
respect to any ISV S m , assuming limited path-dependence [ Schapery (1990)], yields 
the thermodynamic force, / m , available for advancing structural changes associated 
with the m th ISV. 


_ m, 

8S m 


(i) 


It is shown in Schapery (1989, 1990) that the total work potential is stationary with 
respect to each ISV. 


dW T 

~dS^ 


(2) 


Additionally, Rice (1971) utilized the second law of thermodynamics to establish the 
inequality: 

fmSm > 0 (3) 


which suggests that “healing” is not allowed for a material undergoing structural 
changes. Equations (1), (2), and (3) form the foundation of a thermodynamically- 
based work potential theory for modeling non-linear structural changes in a material 
exhibiting limited path-dependence. 
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2.2 Multiple ISV Formulation of ST to Account for Multiple 


Damage and Failure Mechanisms 

Due to the generality of the evolution equations, Equations (2) and (3), the work 
potential theory can account for any number and type of structural changes that 
may occur in a material. This is especially useful for modeling progressive damage 
in composites because the heterogeneity of the composite, and multiaxiality of the 
local fields, enables multiple damage mechanisms to arise during a typical loading 
history. For instance in the matrix phase alone, microdamage accrues until its effects 
are superseded by the growth of larger transverse cracks. Microdamage is considered 
the advancement of microcracks, voids, fissures, shear bands, and other flaws that are 
present in the matrix of a composite [ Sicking (1992); Schapery and Sicking (1995); 
Basu et al. (2006); Ng et al. (2010)]. The size of these flaws is typically on the order 
of that of the fiber or smaller. Transverse cracks nucleate from pre-existing flaws 
within the matrix but grow parallel to the fibers and span the thickness of the lamina 
[ Talreja (1985b); Allen et al. (1987a); Laws and Dvorak (1988); Gudmundson and 
Ostlund (1992); Yang and Cox (2005); Noda et al. (2006); Green et al. (2007)]. Often, 
the growth of individual transverse cracks is extremely rapid; however, the effects of 
transverse cracking on the stiffness of a composite laminate can be progressive if 
multiple cracks form over an extended period of time and throughout an expansive 
volume. Eventually, transverse cracking is succeeded by more catastrophic damage 
mechanisms including interlaminar delamination, fiber breakage, pullout and bridging 
associated with macroscopic laminate fracture [ McCartney (1992a,b); Hallett et al. 
(2008)]. 

The present EST formulation assumes that three major intralaminar mechanisms 
are responsible for all observed non-linearities in the stress-strain curve of a composite 
lamina: matrix microdamage, matrix macroscopic cracking, and axial fiber failure. 
Each of these mechanisms can be accommodated by partitioning the total dissipated 
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energy density, Ws 1 into portions associated with each mechanism. 

Matrix microdamage is the primary cause of observed non-linearity in the stress 
versus strain response of many polymer matrix composites (PMCs) (i.e. systems ex- 
hibiting negligible non-linear elasticity, plasticity or viscous effects) up to localization 
of microdamage into more severe failure mechanisms, such as transverse cracking, 
fiber breakage, kink band formation, or delamination. Microdamage can be consid- 
ered the combination of matrix microcracking, micro-void growth, shear banding, 
and fiber-matrix debonding. Figure 2 shows a typical uniaxial response of a material 
exhibiting microdamage evolution, where the recoverable energy potential is given 
by W and the potential dissipated into evolving structural changes associated with 
microdamage is given by S. 

Typically, matrix microdamage continues to grow until the onset of more catas- 
trophic failure mechanisms initiate. It should be noted that this work explicitly 
distinguishes between damage and failure in the following manner: 

Damage - Structural changes in a material that manifest as pre-peak non-linearity in 
the stress-strain response of the material through the degradation of the secant 
moduli. 

Failure - Structural changes that result from damage localization in a material and 
manifest as post-peak strain softening in the stress-strain response of the ma- 
terial. 

Here, three major failure mechanisms, which are distinct from the microdamage mode, 
are considered: transverse (mode I) matrix cracking, shear (mode II) matrix cracking, 
and axial (mode I) fiber fracture. These failure modes are consistent with the in-plane 
failure typically observed in PMC laminates. It is assumed that the evolution of these 
mechanisms yields an immediate reduction in the load-carrying capability of a local 
subvolume where the mechanism is active. Three IS Vs are used to account for mode 
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I matrix cracking, mode II matrix cracking, and mode I fiber failure, respectively: 
Sy\ and S{. These ISVs are defined completely in Section 2.3, and are taken to 
be the potentials required to advance structural changes associated with these failure 
mechanisms. 

At any given state the total dissipated energy density W$ can be calculated as a 
sum of energy dissipated through the aforementioned damage and failure mechanisms, 
given by the four ISVs. 

Ws = S + S? F + S? IF + s{ (4) 

According to the first law of thermodynamics, the total work potential (ignoring 
thermal dissipation) is given by the sum of the elastic strain energy density and the 
potentials associated with each of the damage or failure mechanisms. 

W T = W E + S + S? + S?j + 5/ (5) 


where IV e is the elastic strain energy density. Invoking the stationarity principle, 
Equation (2), 


dW E 

~dS~ 


dW E 

dS™ 


( 6 ) 


dW E 

dS™ 

dW E 

~ds[ 
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and the Second Law of Thermodynamics, Equation (3), gives: 


S > 0 
ST f > 0 

( 7 ) 

S?IP > 0 

S ! F > 0 

Equations (6) and (7) constitute the evolution equations for damage and failure in 
a material associated with matrix microdamage, matrix cracking, and fiber breakage 
in tension. 

It should be noted, that EST can also account for kink band formation under 
axial compression [Schapery (1995); Basu (2005); Basu et al. (2006)]; although, the 
applied loading in the examples presented in Sections 3 and 4 are tensile, and kink 
banding does not occur. As the lamina is loaded, the fibers in the composite rotate 
by some angle <f>, given by the deformation gradient in the model. To model the kink 
band mechanism, all calculations are then executed in the instantaneous fiber frame 
given by (/>; therefore, fibers rotation induces larger shear strains, 712. Increased shear 
strain yields more damage, leading to a reduction in the shear modulus. The increase 
in shear compliance allows for further progression of the shear strain. Under axial 
compression, this leads to a runaway instability, and a kink band will form. 

2.3 Failure Initiation 

Matrix microdamage requires no initiation criterion. For low strain levels, the micro- 
damage ISV S remains small and its effects on the composite moduli are not apparent. 
As S evolves, with increased strains, its effects on the stress-strain response of the 
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composite become more noticeable. However, it is postulated that the evolution of the 
failure mechanisms immediately yield a negative tangent stiffness; therefore, initia- 
tion criteria are required. Furthermore, criteria are required to mark failure initiation 
because the macroscopic cracks responsible for failure may result from localization 
of microdamage, or they may nucleate from pre-existing flaws in the material not 
necessarily associated with microdamage. 

EST is implemented in homogenized laminae; therefore, phenomenological criteria 
must be utilized that account for the composite microstructure. The Hashin-Rotem 
failure criterion incorporates separate equations for matrix failure and fiber failure 
initiation [ Hashin and Rotem (1973)]. The matrix failure criterion involves contribu- 
tions from both the transverse (€22) and shear (712) strains. 


(|)7f A 

( 8 ) 

(sp (?)’-■ <*<» 


where Y? is the transverse lamina failure strain in tension, Yq is the transverse failure 
lamina strain in compression, and Z is the shear failure strain. The fiber failure 
criterion only involves the axial strain €n. 


£11 

Xj' 


1 €11 > 0 


(9) 


where Xt is the maximium allowable axial strain of the lamina. A local, lamina 
coordinate frame is chosen such that, 1- is the axial direction of the fibers, 2- is the 
in-plane transverse direction, and 3- is the out-of-plane direction. When Equation 

(8) is satisfied the matrix failure ISVs S™ and Sj} are activated, and when Equation 

(9) is satisfied fiber failure evolution Sj is permitted; otherwise, S remains the only 
active ISV. Upon satisfaction of either Equation (8) or Equation (9), it is assumed 
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that the more severe failure mechanisms dominate, superseding the effects of matrix 
microdamage; therefore, S — 0, and additional microdamage is precluded. 

2.4 Use of Traction-Separation Relationships to Define the 
Failure Potentials 

Sicking (1992) and Schapery and Sicking (1995) used a single ISV to model the effects 
of transverse cracking on a composite lamina. Similar to microdamage, the transverse 
and shear moduli were related to transverse crack evolution through a set damage 
functions obtained from coupon experiments. Predictions of the non-linear response 
of numerous laminates were presented assuming a homogenous strain state in the 
laminates. Pineda et al. (2010) implemented the dual-ISV formulation of ST for 
predicting microdamage and transverse cracking within FEM to model the response 
of a center-notched laminate that was reinforced with buffer strips. The original 
formulation required the solution of two, coupled, bi-variate polynomials, which in an 
FEM framework became extremely computationally intensive. Thus, Pineda et al. 
(2010) decoupled the microdamage and transverse cracking evolution equations. 

In the aforementioned publications, it was assumed that the transverse cracking 
affected the relationships between stress and strain. However, the existence of a 
macroscopic crack invalidates the assumption of a continuum. Here, it is presumed 
that failure arises from the evolution of cohesive cracks within the continuum, and the 
ISVs associated with failure (axial, transverse, and shear) influence the relationship 
between traction on the crack faces and the crack-tip opening displacement. The 
satisfaction of Equations (8) and/or (9) indicates the material behavior transitions 
from that of a damaging continuum to that of a cohesive crack, and the essential 
fields become traction and separation, rather than stress and strain (see Figure 3). 

Once a cohesive crack initiates in the continuum, opening of the crack yields a 
reduction in traction on the crack faces at the crack tip. If subsequently the crack 
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is closed, it is assumed that traction at the crack tip will unload linearly towards 
the origin of the traction versus separation law (see Figure 4). The strain energy 
release rate (SERR) G 3 M is taken as the total energy dissipated per unit area of new 
surface that is created through crack advancement and can be calculated as the area 
under the traction-separation law (for a given traction and separation pair) minus 
the energy per area that can potentially be recovered by unloading. 




t J dd 3 


.f) ( \> 

dM u M 


( 10 ) 


where j indicates the material (fiber / or matrix m ) , M represents the corresponding 
mode (mode I or mode II), S J M is the crack tip opening displacement in mode M and 
material j, and t J M is the corresponding traction at the crack tip. 

Theoretically, the shape of the traction-separation laws for mode I crack growth 
in the fiber, and mode I and II crack growth in the matrix can take any shape. 
Gustafson and Waas (2009) investigated triangular, trapezoidal, beta distribution, 
and sinusoidal traction-separation laws in discrete cohesive zone method (DCZM) el- 
ements and determined that the shape only affected convergence of the FEM solver, 
but not the overall results. For simplicity, it is assumed here that all three types of 
cracks obey triangular traction-separation laws, presented in Figure 4. The total area 
under the traction-separation curves is controlled by the corresponding material frac- 
ture toughness in the appropriate mode, where G J IC is the mode I fracture toughness 
of the fiber, Cpc is the mode I fracture toughness of the matrix, and G 7 f l IC is the mode 
II fracture toughness of the matrix. The cohesive strengths of the materials tj C (mode 
I fiber strength), fp c (mode I matrix strength), and tf IC (mode II matrix strength) 
are given by the stresses in the continuum when Equations (8) and/or (9) are satis- 
fied. Mode I, normal cracks are not allowed to grow under compression, but mode II, 
shear cracks can evolve under normal compression. Therefore, the mode I traction- 
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separation laws for the fiber and matrix (Figures 4a and 4b) do not accommodate 
negative crack tip displacements. However under negative mode II displacement (see 
Figure 4c), the traction on the crack faces will increase linearly until the maximum, 
previously attained displacement magnitude is reached, after which, the crack faces 
will resume unloading according to the negative portion of the traction-separation 
law. The traction-separation laws exhibited in Figure 4 do not require any initial, 
fictitious, pre-peak stiffness because the cracks are embedded within a continuum. 
This is an advantage over the use of DCZM elements which do require an initial stiff- 
ness because these interfacial elements do not actually represent physical material 
within the model and must attempt to simulate initially perfect bonding between 
adjacent material domains [Xie et al. (2006); Gustafson (2008)]. If set incorrectly, 
these fictitious stiffnesses can cause numerical problems [ Turon et al. (2006)]. 

Although no mode I crack can advance under compression, it is possible for post- 
peak softening to occur under compressive loading situations. For instance a kink 
band could form under global axial compression, or the matrix could fail in local 
shear due to internal friction (Mohr-Coulomb) in quasi-brittle materials under trans- 
verse compression [Hoek and Bieniawski (1965)]. Since these failure mechanisms 
involve local shear at a the fiber/matrix scales which is typically below the operating 
lamina/laminate scale, it appears that these mechanisms evolve under mode I com- 
pression. In a model containing homogenized laminae, there is no subscale shear to 
drive these compressive failure modes. However, EST could be extended further to 
incorporate these mechanisms through phenomenological accessions. The methods, 
developed by Basu (2005); Basu et al. (2006) and described in Section 2.3, can be used 
to track the instantaneous fiber angle, and a critical fiber angle can be assigned to 
indicate the initiation of post-peak softening due to kink band formation. Similarly, 
a matrix compression criterion, such as the one developed by Buck and Schilrmann 
(1998, 2002) could be used to signal the initiation of Mohr-Coulomb compressive 
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failure. The traction-separation laws for mode I fiber compression and mode I trans- 
verse matrix compression could be adjusted to include the post-peak softening effects 
of microbuckling and Mohr- Coulomb matrix failure. These postulated, compressive, 
mode I traction-separation laws could account for energy released through these sub- 
scale failure modes in a homogenous model at the lamina/laminate scale. However, 
the examples presented in this chapter are tension dominated, and extension of the 
theory to accommodate apparent mode I compressive failure is left for future work. 

Using the traction-separation laws in Figure 4, the SERR can be calculated with 
Equation (10). 

G{ = \t f IC 6{ ( 11 ) 

G? = \tr c 6? (12) 

GTi = \Wic8n ( 13 ) 


It is assumed that the energy released due to cracking is smeared over the entire 
element [ Bazant (1982); Bazant and Oh (1983)]. Thus, the dissipation potentials in 
an element resulting from macroscopic cracking are related to the SERRs using the 
suitable element dimensions. 


nf 

of _ 

I 7(0+90°) 
te 


pm 

nm I 

~ ,00 

he 


£im 1 1 

bn - m 


(14) 

(15) 

(16) 


If there is a single integration point in the element, li 9+90 ^ is the length of a line 
running perpendicular to fiber direction in the element that intersects two edges of 
the element and the integration point, and j is the length of a line that is parallel 
to the fiber direction in the element that intersects two edges of the element and 
the integration point. If there is more than one integration point in the element, 
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the element can be partitioned into a number of subvolumes equal to the number of 
integration points, and the lengths li° +9 ° ^ and le ^ are lengths of lines that intersect 
the corresponding integration point as well as two element edges or integration point 
subvolume boundaries. Incorporating a length scale into the IS Vs results in mesh 
objective, post-peak, softening. This is elaborated upon further in Section 4 . 


2.5 EST Evolution Equations for a Fiber-Reinforced Lamina 

To arrive at the evolution equations for the four ISVs, the elastic strain energy density 
must be defined for a material which may contain cohesive cracks. Therefore, the 
elastic strain energy We is comprised of a contribution from the continuum W and 
any possible cohesive cracks W J M . The plane stress, elastic strain energy density in 
the continuum is defined as 

W — -{Eu^h + E 2 2 {S)e \ 2 + ^12(5)712) + Q 12^11^22 ( 17 ) 

where stress in the laminae are related to strain assuming plane stress conditions. 


<Tli — Qii<ui + Q 12^22 
V22 = Q 12^11 + Q22^22 

Tl2 = Q 66712 


where 712 is the engineering shear strain and 


Q n 
Q22 


En 

1 - ^12^21 
E22 

1 - ^12^21 


Q 12 — V12Q22 


Qm — G12 


^21 


V12E22 

Eu 


( 18 ) 


(19) 
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where Eu is the axial elastic modulus, E 2 2 is the transverse elastic modulus, ^12 is 
the Poisson’s ratio, v 21 is the transverse Poisson’s ratio, and G 12 is the elastic shear 
modulus. After assuming that the quantity ^12^21 << 1 , Equations ( 19 ) simplify, 


Q 11 — En 
Q 22 — E 2 2 


(20) 


Q 12 — ^ 12^22 

Q 66 — Gl2 


Note that only the transverse and shear moduli (E22 and G 12) are functions of S 


since matrix microdamage only accrues in the matrix of the laminae. The Poisson’s 
ratio is assumed to evolve such that the quantity Q12 — E22M2 remains constant; 
however, this restriction can be relaxed if deemed necessary. The degraded moduli 
are related to the virgin moduli (E 22 o and G 120) and the ISV through a set of mi- 
crodamage functions (e s (S) and g s (S)) that are obtained from three uniaxial coupon 
tests [Schapery ( 1989 ); Sicking ( 1992 ); Schapery and Sicking ( 1995 )]. 


Degrading E 2 2 and G 12 exclusively is consistent with the intralaminar damage typi- 
cally observed in PMC laminates. 

The elastic strain energy density of the cohesive cracks are defined as the recov- 
erable energy per unit crack surface area smeared over the entire element. 
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The tractions in Equations (23)-(25) can be related to the secant stiffness’ in the 
traction-separation laws k J M . 

tj — kjSj (26) 


_jjn j^rri^m 


rm 7 m rm 

L II ~~ ^ II °II 


(27) 

(28) 


Hence, the total elastic strain energy density in the continuum is given by 
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Substituting Equation (29) into Equations (6) gives the ISV evolution equations. 
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(33) 


The use of a reduced ISV S r — S 3 has been employed in Equation (30). Sicking 
(1992) has shown that the use of this reduced ISV yields polynomial forms of the 
microdamage functions in Equations (21) and (22). Using the chain rule and the fact 
that 
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by Equations (11)- (16), the cohesive secant stiffnesses are determined. 

H = -J 
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Evaluating the integrals in Equations (37)-(39), while enforcing k° M — 0 when 5 3 M — 
2G j MC results in expressions for k 3 M in terms of S J M . 
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(40) 

(41) 

(42) 


The thermodynamically consistent stiffnesses derived in Equations (40)- (42) can also 
be derived directly from the traction-separation laws using geometry. 

Finally, it is assumed that following failure initiation the strains are related to the 
crack tip opening displacements by 


;< 9+90 °>6n = ;< 9+90O) 4i + S{ 


(43) 
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li e) 712 = lf> 7 i C 2 + 257 , 


(45) 


where ef l5 e§>, and 7 ^ are the strains when Equations (8) and/or (9) are satisfied. 
Equations (43)-(45) imply that the strain in the continuum remains at the values 
obtained when failure initiates, and that any incremental change in the global strain 
after failure initiation is used wholly to advance the crack tip opening displacement. 
To account for changes in the continuum strain after failure initiates, it can be as- 
sumed that the stress state in the cracked body is homogenous and the tractions on 
the crack tip faces are equal to the stresses in the continuum. Then, the strains in 
Equation (29) can be formulated in terms of the cohesive secant stiffnesses and the 
crack tip opening displacement. However, it is assumed that the evolution of strain 
in the continuum is negligible once cohesive cracks form. Equations (43)- (45) can be 
utilized in Equations (40)-(42) to obtain k J M as functions of the global strain at an 
integration point. 
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Once failure initiates, the effects of failure supersede the effects of microdamage 
and evolution of S ceases. The cohesive stiffness in a cracked element is calculated 
using Equations (46)-(48) for a given strain state; then, Equations (26)-(28) and 
(43)- (45) are used to calculate the tractions on the crack tip faces and the crack tip 
opening displacement. It is assumed that the stress state in the integration point 
subvolume of the element is homogenous, and the tractions on the crack tip faces are 
equal to the stresses in the element. Lastly, the axial, transverse, and shear moduli 
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of the element can be calculated [ Bazant and Oh (1983)]: 
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where and (7^ 2 are the degraded transverse and shear moduli, due to microdam- 
age, when Equation (8) is satisfied. 

For visualization purposed in the FEM simulations, degradation parameters are 
defined which relate the current, degraded stiffnesses to their original values upon 
failure initiation. 
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(54) 


The degradation parameter can have a minimum value of zero, which indicates that 
no degradation has occurred, or a maximum value of one, signaling that the corre- 
sponding modulus has been completely diminished. 

The negative tangent stiffness of the stress-strain curve necessary for post-peak 
strain softening to occur imposes a restriction the maximum allowable element size, 
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as shown by Bazant and Oh (1983). 



(55) 


l ^ < min 



(56) 


The analyst must be careful to ensure the dimensions of any failing elements are 
smaller than the conditions given in Equations (55) and (56). 

In summary, Equations (8) and (9) mark the transition from evolving microdam- 
age to failure to macroscopic cracking. Prior to failure initiation, Equation (30) is 


and Sj} remain zero. Equations (21) and (22) are used to calculate the degraded 
transverse and shear moduli. Subsequent to matrix failure initiation, microdamage 
growth is precluded, and S r remains at £*, the value of S r when Equation (8) was 
satisfied. The degeneration of the transverse and shear moduli, resulting from matrix 
transverse and shear cracking, is calculated using Equations (50) and (51). Finally 
if Equation (9) is satisfied, the axial modulus is calculated using Equations (49) as 
fiber breakage evolves in the element. Once the material moduli have been calculated 
using the appropriate evolution equations, the stresses can be updated accordingly 
using Equations (20). 


The theory outlined in Section 2 eliminates the mesh dependency that arises from 
ill-posedness when the elements enter the post-peak softening regime by introducing 
a characteristic length into the formulation. The total SERR dissipated during the 
evolution of the discontinuity is equal to the prescribed fracture toughness and is 
independent of the element size. This approach, commonly referred to as the smeared 


used to calculate the microdamage reduced ISV 5 r , and the failure IS Vs S{, S] 


'm 


3 Mesh Objectivity 
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crack approach, or crack band model, lias been used to alleviate mesh dependency in 
FEM since it was first developed by Bazant (1982) for post-peak strain softening in 
concrete. 

To exhibit the mesh objectivity of EST, a simple example is presented in this 
section. One quarter of a 200 mm x 100 mm panel containing a central hole is 
modeled with finite elements using the Abaqus, version 6.10-1 finite element software 
Abaqus (2008). The panel contains a hole with a radius of — 5 mm in the center. 
The left edge of the panel is constrained in the ^-direction to simulate symmetry. 
Similarly, the bottom edge is constrained in the y-direction. A uniform displacement 
is applied to all the nodes on the top edge of the panel in ^/-direction. Details on 
the panel geometry and boundary conditions are displayed in Figure 5. The panel is 
composed of a generic, [90°], composite lamina with the fiber angle measured with 
respect to the y- axis; thus, the applied displacement is perpendicular to the fiber 
direction in the panel. EST is used to model damage and failure in the panel. 

Four different meshes are used to evaluate the effect of mesh size on the response 
of the panel. All four meshes consist of two-dimensional (2-D), plane stress, quadri- 
lateral, S4R shell elements. The elements are linear, reduced integration elements and 
contain four nodes and one integration point each. The density of the four meshes, 
shown in Figure 6, increases within a region near the central hole. Average element 
sizes equal 0.5^, 0.2?^, O.lr^, and 0.04 are used in the four different meshes. Coarser 
elements are used away from the hole to improve computation time. The four meshes 
are subjected to the same boundary conditions and loading, and are composed of the 
same material properties. The same elastic, damage and failure parameters are also 
used in all four simulations. 

The resultant, applied tensile stress (given by two times the sum of the reaction 
forces at the nodes on the top edge divided by the cross-sectional area) normalized 
by the transverse, mode I, critical strain times the transverse Young’s modulus a is 
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plotted versus the applied displacement normalized by the radius of the hole A for 
the four different meshes in Figure 7. It can be seen that the mesh density has a 
minimal effect on the load-deflection results. The small discrepancy in the results 
between the four meshes can be attributed to the increased accuracy in the fields as 
the mesh is refined. Moreover, the total energy dissipated is preserved from mesh to 
mesh. 

Figure 8 displays contours of the normalized, reduced microdamage ISV S r imme- 
diately before failure initiation in the four different meshes. S r is normalized by the 
maximum S r obtained in all four simulations which is 27% of the maximum allowable 
S r required to bring the moduli to zero. The four meshes exhibit similar microdamage 
contours, but as the mesh is refined, the microdamage is contained to the vicinity of 
the hole. Additionally, the global stress at which failure initiates reduces slightly as 
the mesh is refined; this supports the previous statement that the discrepancies in 
the stress-displacement responses were a facet of increasing held accuracy with mesh 
refinement. 

The transverse matrix failure degradation parameter DJ 1 is plotted for all four 
meshes in Figure 9. Figures 9a-9d show the failure pattern at the ultimate, global 
load. The coarsest mesh shows that a crack band has grown at the intersection of 
the hole and the bottom symmetric boundary and is propagating towards the free 
edge, while moving away from the bottom boundary when the specimen ultimate 
load is achieved. In the hirer meshes, the crack band path is different. In Figures 
9b-9c the crack band initiates at the hole slightly above the bottom boundary and 
propagates towards the free edge while approaching the bottom boundary. The hnest 
mesh, Figure 9d, exhibits multiple crack bands near the hole, but only one of the 
crack bands grows signihcantly. While in all simulations the hrst crack band initiates 
at the hole, Figures 9c and 9d show some crack bands beginning to initiate at the free 
edge and propagate inwards by the time the ultimate load is attained. Figures 9e-9h 
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display the transverse crack band pattern once the specimen has lost all load-carrying 
capability. The solution for the simulation with the finest mesh 0.04?^ diverged before 
all the load carrying capability was lost; so, Figure 9h presents the crack band path 
at the final converged state, which is still far below the ultimate load state. In 
Figures 9e-9h the same crack band patterns that developed in Figures 9a-9d are 
evident, except those crack bands have saturated to maximum degradation: D™ — 1. 
Additionally, the crack bands advancing from the free edge observed in Figures 9c and 
9d have progressed further. However, this is well after the specimens reached their 
ultimate load; therefore, it is assumed that those crack bands did not influence the 
load carrying capability of the panels. The discrepancy in crack band path observed 
for the different meshes indicate that the crack band path is dictated by the accuracy 
of the fields surrounding the leading boundary of the crack band path. 

4 Example - Center Notched Panels Subjected to 
Uniaxial Tension 

4.1 Experimental Details 

Two center-notched panel (CNP) configurations were tested at the NASA Langley 
Research Center (LaRC) [Bogert et al. (2006); Satyanarayana et al. (2007)]. The 
geometrical details of the panel and testing boundary conditions are presented in 
Figure 10. The panels were 3” wide and 11.5” long. Two 3” x 2.75” tabs were placed 
on both ends of the specimens, leaving a gage section of 3” x 6” which is displayed 
in Figure 10. A central notch was machined in each panel that was 0.75” wide and 
had a notch tip radius of 0.09375”. The end tabs were clamped and a vertical, 
tensile displacement (in the y-direction) was applied to the top tab using a servo- 
hydraulic testing machine. The bottom tab was fixed preventing any ^-displacement 
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of the bottom boundary of the gage section. The gripped tabs also prevented any 
displacement in the ^-direction at the top and bottom boundaries of the gage section. 

The panels were comprised of laminated T800/3900-2 carbon fiber/toughened 
epoxy composites. Three different lay-up configurations were tested; however one of 
the configurations exhibited significant delamination. Since the focus of this work 
is modeling in-plane damage and failure mechanisms, this configuration is not con- 
sidered here. The two remaining configurations are presented in Table 1. The first 
lay-up, Laminate-1, consists of 12 0° plies, and the second, Laminate-2, is a symmet- 
ric, multi-angle lay-up with 40% |45°|, 40% 0°, and 20% 90° layers. 

Several strain gages where affixed to the test panel, labeled Sg-1 through Sg-4 in 
Figure 10. Sg-1 was placed in the center of the panel, 1.5” above the notch. Sg-2 was 
placed 1.5” above the notch tip. Sg-3 was attached in front of the notch, 0.5” from the 
free edge, and Sg-4 placed at the notch tip. Global load versus displacement data, and 
local strain gage data was reported in Refs. [Bogert et al. (2006); Satyanarayana et al. 
(2007)], along with a post-test C-Scan of Laminate-1 and photograph of Laminate-1 
and photograph of Laminate- 2. 

4.2 Finite Element Model Details 

The linear elastic properties of T800/3900-2 used in the FEM models are presented 
in Table 2, and were taken from Ref. [Bogert et al. (2006)]. The shear microdam- 
age function g s utilized in Equation (22) was obtained from [45 0 /-45°]3s< angle-ply 
T800/3900-2 coupon tests. The transverse, tensile and compressive microdamage 
functions were inferred by scaling the coefficients of the microdamage curves pre- 
sented by Sicking (1992) for AS4/3502 by the ratio of the virgin transverse modulus 
of T800/3900-2 to that of AS4/3502, as the stress-strain curves of the coupon lami- 
nates necessary to characterize e s were not available. The polynomial forms of e s and 
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g s are 


G s (S r ) — e S Q + e s i S r + e S 2 + e s s + e S 4 + e s ^S^ 


(57) 


g 8 (S r ) — g s 0 + 9 s \Sr + 9s2$r + 9s3^r + + 0s5$r 


(58) 


The coefficients of the microdamage curves are presented in Table 3, and the curves 
are plotted in Figure 12. 

To increase computational efficiency, the first derivatives of the higher order mi- 
crodamage polynomials are approximated using linear spline interpolants. Thus, the 
microdamage evolution equation, Equation (30), is always second order in S r yielding 
a very efficient analytical solution. Since the value of S r from the previous increment 
is used to estimate which spline regime should be used to solve for S r in the current 
increment, the solution is checked to ensure that it falls within the applicable range 
of the spline that was used. If it falls outside of the range of S r that are valid for 
the splines, the solution is calculated again using splines that accord to the solution 
of the previous iteration. This procedure continues until the solution of Equation 
(30) falls within the relevant range of S r for the splines used in Equation (30). A 
maximum number of iterations can be set, after which Equation (30) is solved using 
the full polynomial forms of the damage functions, given in Equations (57) and (58), 
by finding the eigenvalues of the companion matrix of the polynomial coefficients of 
the evolution equation. 

The axial mode I, transverse mode I, and shear mode II critical cohesive strains, 
and fracture toughness’ are given in Table 4. The matrix mode I and mode II cohesive 
critical strains (Yt, Yc , and Z) and the fracture toughnesses ( G™ c and G™ IC ) were 
calibrated using data from Laminate- 1. These values were adjusted until the global 
load versus local strain at strain gage Sg-1 (see Figure 10) obtained from the model 
provided the best match to the experimental data. Both the simulation results and 
experimental data are presented in Figure 14a. Laminate- 1 did not exhibit any axial 
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failure; so, the fiber mode I parameters (Xt and G{ c ) were calibrated such that 
the ultimate load from the simulation of Laminate-2 corresponded with the ultimate 
load reported by Bogert et al. (2006) for Laminate-2. Subsequent work will outline a 
procedure for measuring these values experimentally. 

Displacement was applied to both laminates using the *DYNAMIC keyword in 
Abaqus with the parameter APPLICATION = QUASI-STATIC. This implicit dynamic 
solver is recommended for quasi-static problems exhibiting a high-degree of nonlinear- 
ity. This procedure uses numerical damping to stabilize the problem. The numerical 
damping does not significantly affect the simulation results because the velocities in 
these simulations are low. For Laminate-1 a total displacement of 0.0236”, and for 
Laminate-2 a displacement of 0.0472”, is applied over 1000 seconds. The panels were 
assigned a representative density of 0.057 lb/in. 3 . This implicit, dynamic technique 
has advantages over traditional static, implicit solvers which have difficulty converging 
when the material exhibits post-peak softening [Belytschko et al. (2000); Belytschko 
and Mish (2001)], and is not limited by a minimum stable time step required with 
explicit solvers [ Hughes (2000)]. 

4.3 Results - Laminate-1 

Global load P versus displacement A of a 4” section of Laminate- 1 is compared 
to results from the EST simulation in Figure 13. Very good agreement between the 
model and the experimental results is achieved. The response of the specimen appears 
to be linear until near 8,000 lbf. , where the specimen begins deforming non- linearly. 
The EST simulation captures the initiation and progression of the global nonlinearity 
accurately. This panel was not loaded until catastrophic failure; hence, the data 
presented in Figure 13 represents load versus displacement data prior to the ultimate 
load of the specimen. 

Local strain gage data (global load P versus local ^-direction strain e yy ) from 
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Laminate-1 is plotted with the results from the EST FEM model in Figure 14; please 
refer to Figure 10 for locations of strain gages. Strain relaxation is observed in the gage 
farthest away from the notch: Sg-1, shown in Figure 14a. The mode I and mode II 
matrix failure parameters in EST were calibrated such that the model demonstrates 
the same transition into strain relaxation at this location and at a similar global 
applied load. This load, taken as the splitting load, is 8,250 lbf. in experiment and 
8,210 lbf. in the model (summarized in Table 5). The transition to strain relaxation 
is more abrupt in the experiment as evidenced by the sharp knee in the load-strain 
curve, whereas, the transition in the model is more gradual. Prior to strain relaxation 
at this point, the experiment displayed slight stiffening not observed in the model. 
Additionally, the model response is much smoother than that of the experiment in 
the strain relaxation regime. Even though the global loading is quasi-static, local 
events, such as cracking, may be dynamic; therefore, the discrepancy in the strain 
relaxation portion of the load-strain curves could be a result of dynamic matrix crack 
growth and arrest in the test specimen. Local crack dynamics were not taken into 
account in the model. Additionally, the jaggedness of the experimental data may be 
a facet of the stochastics related to the local microstructure of the composite that 
are not included in the model. The data from the experiment and simulation for 
Sg-2, which is located 1.5” directly above the notch, are presented in Figure 14b. 
The model predicts less strain at Sg-2, for a given load, than the experiment, but the 
non-linear trends are very similar. This gage lies directly in front of the splitting crack 
path, shown in Figure 10, and it is not realistic to expect perfect agreement in areas 
experiencing high levels of damage and failure because of idealizations used to model 
the evolution of cracks in the simulations. Figure 14c displays data for Sg-3, located in 
front of the notch near the free edge. Very good agreement between the experimental 
and simulation results are exhibited. The model accurately captures the non-linear 
evolution of strain, away from the highly damaged regions, as a function of applied 
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load. Finally, results for Sg-4 (located directly at the notch tip) are given in Figure 
14d and includes the experiment and simulation display of axial strain relaxation. 
As with Sg-2, Sg-4 shows less strain for a given applied load. However, the load at 
which the strain at Sg-4 relaxes in both the experiment and model correlate well, in 
accordance with the splitting load. Again, the relaxation response of the experiment 
is discontinuous, but the model exhibits continuous behavior. 

A C-Scan of the failed Laminate- 1 specimen is displayed in Figure 15. Four 
splitting cracks can be observed propagating outward from the notch tip, parallel to 
the loading direction, towards the gripped edges. Contour plots of the normalized 
microdamage obtained from the simulation are presented in Figure 16 at the splitting 
load 8,200 lbf. and at 16,400 lbf. In these plots, S r is normalized by the maximum 
achievable value, S r r nax — 7.57 psi^, obtained from Figure 12. In Laminate-1, S r 
reached a maximum value equal to 0.171 S™ ax . At the splitting load, the regions of 
maximum damage are localized to small regions, along the same crack path observed 
in Figure 15, embedded in a more widespread domain and exhibiting less severe 
microdamage. A similar microdamage contour is observed at P — 16,400 lbf., except 
the localized damage region has nearly proceeded to the fixed boundaries of the panel. 

Figure 17 shows the shear failure degradation factor DJ} at the splitting load and 
16,400 lbf. The shear failure localizes into crack bands that are a single element wide 
and progress equivalently to the cracks observed in the experiment. At the splitting 
load (Figure 17a), the crack bands have progressed less than half of the way between 
the notch and the panel boundary on either side of the notch. In Figure 17b, the 
crack bands have nearly reached the fixed grip boundaries. Additionally, Figure 17b 
displays some irregularity in the crack path. In these regions, the mesh is not uniform 
because it needed to accommodate larger elements used to represent the strain gages 
(see Figure 11). No axial failure was observed in this simulation. The shear crack path 
was in a state of transverse compressions, and according to the traction-separation 
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relations used in Figure 4b, transverse failure does not progress under compressive 
conditions. Thus, contours of D j and D ™ are not shown. 

4.4 Results - Laminate-2 

Numerical results for applied load versus displacement of a 4” section of Laminate-2 
are presented in Figure 18. The experimental ultimate load 15,300 lbf. correlates well 
(axial failure parameters were calibrated to obtain an ultimate load that most closely 
matched the experimental data) with the ultimate load obtained from the model, also 
15,300 lbf., and is summarized with the splitting load from the Laminate-1 analysis in 
Table 5. The global response up to failure is nearly linear and failure occurs suddenly 
and catastrophically. 

Figure 19 compares the applied load versus strain gage results from the model 
to the data from the experiment. Sg-1 and Sg-2 exhibited similar behavior; the 
strain increases until the ultimate load is obtained, after which the strain relaxes 
abruptly. The experimental data and numerical results both display this behavior. 
The model exhibits slightly more strain, for a given load, prior to ultimate failure. 
At Sg-3, the model predicts strain localization after the ultimate load is achieved. 
The gage data shows a slight reduction in strain as the load drops; however, the gage 
was placed directly in the crack path and may have been damaged when the panel 
failed. The model results and experimental data for Sg-4 exhibit similar trends, but 
the strain gage shows a large degree of nonlinearity at the notch tip. Bogert et al. 
(2006) attributed this observed nonlinearity to local interlaminar stresses near the 
notch free edge which caused some local delaminations. Since the focus of this work 
was modeling in-plane damage mechanisms, these effects are not captured; however, 
the model could be easily extended to incorporate delamination by placing DCZM 
elements between continuum shell layers [ Satyanarayana et al. (2007)]. 

In the experiment, the gages measure the strain over a continuous area associated 


NAS A/TM— 20 11-217401 


30 


with the size of the gage, but in the model, the strain is taken at the integration point 
of an element; thus, these measures should not be expected to correspond exactly. 
In areas where there are large gradients present, such as near a notch tip (Sg-4) or 
near cracks (Laminate- 1, Sg-2 or Laminate-2, Sg-3), it becomes even more difficult to 
relate the strain gage data to numerical strains from a discretized continuum. This 
may contribute to some of the discrepancies between the local gage data and the 
model results in Figures 14 and 19. 

A photograph taken of the failed, Laminate- 2 specimen is presented in Figure 20. 
The photograph shows that two macroscopic cracks initially propagate from the notch 
tip towards the free edges, perpendicular to the applied load, in a self-similar fashion. 
Eventually, the cracks turn and proceed towards the free edge at an angle. Bogert 
et al. (2006) claim, supported by visual image correlation displacement data, that 
there was some eccentricity in the specimen alignment, which resulted in deviation 
from self-similar crack growth. 

Normalized microdamage contours just prior to the ultimate load are presented for 
the outermost 45°, 0°, -45°, and 90° plies in Figure 21. Similar microdamage patterns 
are evident in the 45° and -45° layers. Microdamage propagates outward, toward the 
free edge, from the notch tip in petal-like patterns. The microdamage in these layers 
is highly distributed throughout the plies. The 0° ply displays a more contained 
microdamage pattern along the lines of the microdamage contours associated with 
axial splitting as shown in Figure 16. A moderate level of microdamage is also 
displayed in the 90° layer, but a low degree of microdamage is distributed throughout 
most of the layer. 

Figure 22 shows the axial failure degradation parameter Dj at the ultimate load for 
the four unique layers. A small amount of axial failure in the 45°, 0°, and -45° layers 
can be observed at the notch tips. It appears that more failure occurs at one notch 
tip than the other. This can be attributed to numerical imperfections resulting from 
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dissimilar meshes at the opposite notch tips, that is the mesh is not symmetric about 
the y- axis. No axial failure is observed in the 90° layer. Contours of the transverse, 
mode I, failure degradation parameter D at the ultimate load are plotted in Figure 
23. The failure patterns are similar in the 45° and -45° plies in Figures 23a and 23c 
and are comparable to the microdamage contours in Figures 21a and 21c, except 
the failure is restricted to regions on either side of the notch. Furthermore, small, 
highly degraded domains can be observed propagating from the notch tip at an angle 
corresponding to the fiber direction in the ply. The 90° layer exhibits some moderate 
degradation in a localized region around the notch tips, and the 0° layer does not 
exhibit much D™. Contours of the shear, mode II, failure degradation parameter D 
are presented at the ultimate load in Figure 24. Very similar failure paths can be 
seen in the 45° and -45° layers and the patterns are nearly symmetric across both 
centerlines of the panel. This is expected because as Figure 4c indicates, the sign of 
the local shear strain does not affect the failure degradation. D in the 0° and 90° 
is limited to very small regions surrounding the notch tips. 

Contours representing the microdamage in the four unique layers are presented 
in Figure 25 after the panel has completely failed and lost all of its load carrying 
capability. Although further matrix microdamage evolution is prohibited in elements 
that have failed (transverse/shear or axial), in the other elements that have not failed, 
matrix microdamage evolution continues. Nearly the entire 45° and -45° layers reach 
a microdamage level of 0.18 S™ ax . The 0° and 90° plies exhibit similar microdamage 
patterns; however, low levels of microdamage are more widespread in the 90° ply. 
Figure 26 shows the fiber failure path once the specimen has completely failed. All of 
the layers, except the 90° layer, show self similar cracks propagating from the notch 
tips towards the free edges of the panel. The angled crack path shown in Figure 
20 was not reproduced because the eccentric loading (suspected in the test) was not 
introduced into the simulation; therefore, the crack growth remained self-similar. A 


NAS A/TM— 20 11-217401 


32 


high degree of transverse matrix failure can be seen in the axial crack path in the 
45°, -45°, and 90° plies in Figure 27. In the 0° layer, some transverse failure is 
observed surrounding the fiber failure, as well as away from the axial failure path, 
which resulted from a stress wave reflecting off the free edges when the axial crack 
band reaches the boundary. Finally, D y} is presented after the specimen has failed 
in Figure 28. Similar failure to Figures 24a and 24c in the 45° and -45° is exhibited, 
but a highly degraded region has localized in the axial crack path. Figures 28b and 
28d show fairly extensive regions containing a high degree of shear matrix failure 
surrounding the axial failure path. 

5 Conclusions 

A thermodynamically-based, work potential theory for damage and failure in compos- 
ite materials, enhanced Schapery theory (EST), was developed. A marked distinction 
between damage and failure was introduced. Damage was considered to be the evo- 
lution of mechanisms that cause structural changes in the material such that the 
non-linear tangent stiffness tensor remains positive definite. Failure was taken to be 
the effect of structural changes in the material that result in loss of positive definite- 
ness of the tangent stiffness matrix and post-peak strain softening. Separate internal 
state variables (IS Vs) were used to account for damage and three in-plane failure 
mechanisms. 

In EST, matrix microdamage, which includes matrix microcracking, shear band- 
ing, and microvoid growth, is responsible for all damage in a composite lamina and 
was accounted for with a single ISV, along the lines of the original Schapery the- 
ory (ST) formulation. The relationship between the transverse and shear moduli 
of the lamina were related to the ISV through a pair of experimentally-obtainable 
microdamage functions. 
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Three major, in-plane failure mechanisms applicable to continuous fiber-reinforced, 
laminated, polymer matrix composites were identified: mode I matrix cracks, mode 
II matrix cracks, and fiber breakage. A failure initiation criterion was used to mark 
the transition from a damaging continuum to a damaged continuum with an embed- 
ded discontinuity. After failure initiation, microdamage evolution ceases and sepa- 
rate IS Vs are introduced to incorporate the effects of the three major failure mech- 
anisms. Evolution of the failure IS Vs is based upon traction-separation laws (which 
are a functions of the appropriate fracture toughnesses) and a characteristic element 
length. Typically, the existence of a non-positive definite stiffness tensor would re- 
sult in pathologically mesh dependent solutions; however, in EST, mesh objectivity 
is ensured by incorporating a characteristic length scale into the failure evolution. 

In Section 3 of this paper, mesh objectivity is demonstrated. A unidirectional 
composite plate with a central hole obeying EST is loaded in transverse tension and 
the response is calculated using four different, increasingly refined, meshes. The global 
stress versus strain response remained unaffected by the change in mesh, save for the 
effects from increased accuracy of local fields in the vicinity of the hole with denser 
meshing. 

Two center-notched panels, with different lay-ups, composed of T800/3900-2 were 
tested under tensile loading at NASA LaRC. Global load versus displacement and 
global load versus local strain gage strain data were compared to results obtained 
from FEM models utilizing EST in Section 4. Quantitatively, very good correlation 
was achieved for both laminates. Furthermore, damage and failure paths predicted 
by the models matched well with the experimental results. 
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ID 

Stacking Sequence 

Thickness (in.) 

Laminate- 1 

[0°] 12 

0.078 

Laminate-2 

[45°/0 o /-45 o /0°/90 o ]s 

0.065 


Table 1: T800/3900-2 lay-up configurations used in CNP tests at NASA LaRC. 


Property 

T800/3900-2 

Ei i (Msi) 

23.2 

E 22 (Msi) 

1.3 

G \2 (Msi) 

0.9 

M2 

0.28 


Table 2: Linear elastic properties for T800/3900-2 used in FEM models. 


e s (S r ) Coefficients 

Tension 

Compression 

g s (S r ) Coefficients 



1.0000 

1.0000 

9so 

1.0000 

1 

-6.0373E-2 

8.4887E-4 

9s l 

-3.0567E-2 

e S 2 

2.5937E-2 

2.8002E-2 

9s 2 

-1.2135E-1 


-1.5789E-2 

-6.2122E-3 

9s 3 

3.7438E-2 

&s4 

2.2571E-3 

N/A 

9s4 

-4.5405E-4 

e S 5 

-1.0440E-4 

N/A 

9s5 

1.9532E-4 


Table 3: Microdamage function coefficients for T800/3900-2 used in FEM models. 


Property 

T800/3900-2 

Xj' 

0.021 

Y t 

0.0048 

Y c 

0.0119 

Z 

0.0075 

Wc 

1026 lbf T' 

m. 

U IC 

170.0 ‘“i 11 - 

m. 

Mm 

^ IIC 

13.54 ‘“i 11 - 

m. 


Table 4: Failure parameters for T800/3900-2. 



Type 

Experimental 

Numerical 

Laminate- 1 

Splitting 

8,250 lbf. 

8,210 lbf. 

Laminate- 2 

Ultimate 

15,300 lbf. 

15,300 lbf. 


Table 5: Critical experimental and simulation loads for Laminate- 1 and Laminate-2. 
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Figure 1: Typical stress-strain curve, containing pre-peak nonlinearity and post-peak 
strain softening, showing the total elastic (We) and total dissipated (Ws) potentials. 



Figure 2: Typical stress-strain curve with a positive-definite tangent stiffness exhibit- 
ing microdamage, showing the elastic (W) and irrecoverable (S) portions. 
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8 22 ’ °22 



Figure 3: Schematic showing the transition form a continuum to a cohesive zone due 
to the initiation of macroscopic cracks. The essential, constitutive variables switch 
from stress and strain to traction and separation. 




(b) Mode I matrix fracture. 



(c) Mode II matrix fracture. 


Figure 4: Triangular traction versus separation which dictates the behavior of cohe- 
sive cracks embedded in the continuum. The total area under the traction-separation 
law represents the material fracture toughness G 3 mC . The area above the unloading 
line for a given traction-separation state is the strain energy release rate G 3 m . 
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Figure 5: Example problem used to demonstrate mesh objectivity of EST. One 

quarter of a 200 mm x 100 mm containing a central hole with a radius of 5 mm 
is loaded in tension with symmetric boundary conditions on the bottom and left 
boundaries. 



(b) 0.2rv 




(c) 0.1 r h . 


(d) 0.047V 


Figure 6: Four mesh densities used to demonstrate the mesh objectivity of EST. 
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Figure 7: Reaction stress normalized by critical axial strain times axial Young’s 

modulus verse applied displacement normalized by hole radius for four different mesh 
densities. 


Reld-2 

SNEG, (fraction = 
(Avg: 75%) 


4-1.000O+00 
4-9. 1670-01 
4-8.333e*01 
4-7. 500e-01 
4-6.6670-01 
4-5.8330-01 
4-S.OOOo-Ol 
4-4. 1670-01 
4-3.3330-01 
4-2.500O-01 


- 1 . 0 ) 





(a) 0.5 Th, cr = 0.50. 


(b) 0.2 r h , a = 

0.42. 


(c) 0.1 r h , a = 

0.37. 


(d) 0.04^, a = 
0.36. 


Figure 8: Contours of the reduced microdamage ISV S r , normalized by the maxi- 

mum S r obtained from all simulations, immediately prior to failure initiation for four 
different mesh densities. 
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Fleld-2 

SNEG, (fraction = 
(Avg: 75%) 

+1.000e+00 
+9.167e-01 
+8.333e-01 
+7.5006-01 
+6.667e-01 
+5.833e-01 
+5.0006-01 
+4. 167e-01 
+3.333e-01 
+2.5006-01 
+1.667e-01 
+8.333e-02 
+0.000e+00 


1.0) 


Y 



X 



(a) 0.5r/j,, a — 0.80. 



(c) 0.1 rhi (7 — 0.82. 



(e) 0.5 r h , a = 0.030. 



(g) 0.1 r h , a = 0.042. 



(b) 0.2r)j, a — 0.81. 



(d) 0.0477,, a = 0.82. 



(f) 0.277,, a = 0.018. 



(h) 0.0477,, a = 0.242. 


Figure 9: Contours of the transverse degradation parameter D ™ , indicative of the 

transverse crack path in the specimens. The contours (a-d) are presented at the 
ultimate load achieved by the specimens and (e-h) after the specimens have lost their 
load carrying capability. 
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Figure 10: Geometry, boundary conditions, and strain gage (Sg) locations of CNPs 
tested at NASA LaRC [Bogert et al. (2006)]. 



Figure 11: FEM mesh used to simulate tensile loading of CNPs. 
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(a) Shear microdamage function obtained 
from [±45°]s angle-ply laminate. 


Figure 12: Microdamage functions 



(b) Transverse tension and compression 
nricrodamage functions obtained by scal- 
ing data for AS4/3502 in Ref. [ Sicking 
(1992)]. 


T800/3900-2 used in FEM models. 



Figure 13: Applied load versus displacement of a 4” section for Laminate- 1. 
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(a) Sg-1. (b) Sg-2. 




(c) Sg-3. (d) Sg-4. 


Figure 14: Applied load versus local strain for Laminate-1. 


Notch 


\ 



Damage Path (Splitting) 


Figure 15: C-Scan of failed Laminate-1 specimen [Bogert et al. (2006)]. 


NAS A/TM— 20 11-217401 


43 


Field- 1 


SNEG, (fraction = -1.0) 
(Avg: 75%) 


+1.714e-01 
+1.571e-01 
+1.428e-01 
+1.285e-01 
+1. 143e-01 
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+8. 569e-02 
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+1.428e-02 
+0.000e+00 




(a) P = 8,210 lbf. 
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Field-2 


SNEG, (fraction = -1.0) 
(Avg: 75%) 


+l,000e+00 
+9. 167e-01 
+8.333e-01 
+7. 500e-01 
+6.667e-01 
+5.833e-01 
+5.000e-01 
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(a) P = 8,210 lbf. 



(b) P = 16,400 lbf. 

Figure 17: Matrix shear failure degradation D in Laminate- 1 
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Figure 18: Applied load versus displacement of a 4” section for Laminate-2. 





(b) Sg-2. 



(d) Sg-4. 

Figure 19: Applied load versus local strain for Laminate-2. 
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Figure 20: Photograph of failed Laminate-2 specimen [Bogert et al. (2006)]. 


Field-1 


LAYER_1 (middle) 
(Avg: 75%) 


i 


+3.147e-01 
+2.884e-01 
+2.622e-01 
+2.360e-01 
+2.098e-01 
+ 1.835e-01 
+1.573e-01 
+ 1.311e-01 
+ 1.049e-01 
+7.866e-02 
+S.244e-02 
+2.622e-02 
+0.000e+00 




(a) 45° Layer. 



(b) 0° Layer. 



(c) -45° Layer. (d) 90° Layer. 


Figure 21: Normalized matrix microdamage contour s ^ ax in Laminate-2 just prior 

to first axial failure initiation P — 8,640 lbf. 
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Field-2 


SNEG, (fraction = -1.0) 
(Avg: 75%) 
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+6.667e-01 
+5.833e-01 
+5.000e-01 
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(a) 45° Layer. (b) 0° Layer. 



(c) -45° Layer. (d) 90° Layer. 


Figure 22: Fiber failure degradation Dj in Laminate-2 at ultimate load P — 15,300 
lbf (magnified view of region near notch). 
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Field-2 


SNEG, (fraction = -1.0) 
(Avg: 75%) 
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(a) 45° Layer. 



(b) 0° Layer. 



(c) -45° Layer. (d) 90° Layer. 

Figure 23: Transverse matrix failure degradation D™ in Laminate-2 at ultimate load 
P = 15,300 lbf. 
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Field-2 
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(a) 45° Layer. 



(b) 0° Layer. 



(c) -45° Layer. (d) 90° Layer. 


Figure 24: Shear matrix failure degradation D y} in Laminate-2 at ultimate load P 

= 15,300 lbf. 
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(a) 45° Layer. 



(b) 0° Layer. 



(c) -45° Layer. (d) 90° Layer. 

Figure 25: Normalized matrix microdamage contour s ^ ax in Laminate-2 after spec- 
imen has lost load carrying capability. 
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(a) 45° Layer. 



(b) 0° Layer. 



(c) -45° Layer. (d) 90° Layer. 


Figure 26: Fiber failure degradation Dj in Laminate- 2 after specimen has lost load 
carrying capability. 
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(a) 45° Layer. 


(b) 0° Layer. 



Figure 27: Transverse matrix failure degradation D™ in Laminate-2 after specimen 
has lost load carrying capability. 
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(a) 45° Layer. 


(b) 0° Layer. 



(c) -45° Layer. (d) 90° Layer. 


Figure 28: Shear matrix failure degradation D™j in Laminate-2 after specimen has 

lost load carrying capability. 
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